{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## MatrixTable Tutorial\n",
    "\n",
    "If you've gotten this far, you're probably thinking:\n",
    "\n",
    "- \"Can't I do all of this in `pandas` or `R`?\" \n",
    "- \"What does this have to do with biology?\"\n",
    "\n",
    "The two crucial features that Hail adds are _scalability_ and the _domain-specific primitives_ needed to work easily with biological data. Fear not! You've learned most of the basic concepts of Hail and now are ready for the bit that makes it possible to represent and compute on genetic matrices: the [MatrixTable](https://hail.is/docs/0.2/hail.MatrixTable.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the last example of the [Table Joins Tutorial](https://hail.is/docs/0.2/tutorials/06-joins.html), the ratings table had a compound key: `movie_id` and `user_id`.  The ratings were secretly a movie-by-user matrix!\n",
    "\n",
    "However, since this matrix is very sparse, it is reasonably represented in a so-called \"coordinate form\" `Table`, where each row of the table is an entry of the sparse matrix. For large and dense matrices (like sequencing data), the per-row overhead of coordinate reresentations is untenable. That's why we built `MatrixTable`, a 2-dimensional generalization of `Table`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### MatrixTable Anatomy\n",
    "\n",
    "Recall that `Table` has two kinds of fields:\n",
    "\n",
    "- global fields\n",
    "- row fields\n",
    "\n",
    "`MatrixTable` has four kinds of fields:\n",
    "\n",
    "- global fields\n",
    "- row fields\n",
    "- column fields\n",
    "- entry fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Row fields are fields that are stored once per row. These can contain information about the rows, or summary data calculated per row.\n",
    "\n",
    "Column fields are stored once per column. These can contain information about the columns, or summary data calculated per column.\n",
    "\n",
    "Entry fields are the piece that makes this structure a matrix -- there is an entry for each (row, column) pair."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Importing and Reading\n",
    "\n",
    "Like tables, matrix tables can be [imported](https://hail.is/docs/0.2/methods/impex.html) from a variety of formats: VCF, (B)GEN, PLINK, TSV, etc.  Matrix tables can also be *read* from a \"native\" matrix table format.  Let's read a sample of prepared [1KG](https://en.wikipedia.org/wiki/1000_Genomes_Project) data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from bokeh.io import output_notebook, show\n",
    "\n",
    "import hail as hl\n",
    "\n",
    "output_notebook()\n",
    "\n",
    "hl.utils.get_1kg('data/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mt = hl.read_matrix_table('data/1kg.mt')\n",
    "mt.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "There are a few things to note:\n",
    "\n",
    "- There is a single column field `s`. This is the sample ID from the VCF. It is also the column key.\n",
    "- There is a compound row key: `locus` and `alleles`.  \n",
    "  - `locus` has type `locus<GRCh37>`\n",
    "  - `alleles` has type `array<str>`\n",
    "- GT has type `call`.  That's a genotype call!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Whereas table expressions could be indexed by nothing or indexed by rows, matrix table expression have four options: nothing, indexed by row, indexed by column, or indexed by row and column (the entries).  Let's see some examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mt.s.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mt.GT.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### MatrixTable operations\n",
    "We belabored the operations on tables because they all have natural analogs (sometimes several) on matrix tables.  For example:\n",
    "\n",
    " - `count` => `count_{rows, cols}` (and `count` which returns both)\n",
    " - `filter` => `filter_{rows, cols, entries}`\n",
    " - `annotate` => `annotate_{rows, cols, entries}` (and globals for both)\n",
    " - `select` => `select_{rows, cols, entries}` (and globals for both)\n",
    " - `transmute` => `transmute_{rows, cols, entries}` (and globals for both)\n",
    " - `group_by` => `group_{rows, cols}_by`\n",
    " - `explode` => `expode_{rows, cols}`\n",
    " - `aggregate` => `aggregate_{rows, cols, entries}`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Some operations are unique to `MatrixTable`:\n",
    "\n",
    "- The row fields can be accessed as a `Table` with [rows](https://hail.is/docs/0.2/hail.MatrixTable.html#hail.MatrixTable.rows)\n",
    "- The column fields can be accessed as a `Table` with [cols](https://hail.is/docs/0.2/hail.MatrixTable.html#hail.MatrixTable.cols).\n",
    "- The entire field space of a `MatrixTable` can be accessed as a coordinate-form `Table` with [entries](https://hail.is/docs/0.2/hail.MatrixTable.html#hail.MatrixTable.entries). Be careful with this! While it's fast to aggregate or query, trying to write this `Table` to disk could produce files _thousands of times larger_ than the corresponding `MatrixTable`.\n",
    "\n",
    "Let's explore `mt` using these tools.  Let's get the size of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mt.count()  # (rows, cols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's look at the first few row keys (variants) and column keys (sample IDs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mt.rows().select().show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mt.s.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's investigate the genotypes and the call rate.  Let's look at the first few genotypes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mt.GT.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's look at the distribution of genotype calls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mt.aggregate_entries(hl.agg.counter(mt.GT.n_alt_alleles()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's compute the overall call rate directly, and then plot the distribution of call rate per variant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "mt.aggregate_entries(hl.agg.fraction(hl.is_defined(mt.GT)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Here's a nice trick: you can use an aggregator inside `annotate_rows` and it will aggregate over columns, that is, summarize the values in the row using the aggregator.  Let's compute and plot call rate per variant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mt2 = mt.annotate_rows(call_rate=hl.agg.fraction(hl.is_defined(mt.GT)))\n",
    "mt2.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "p = hl.plot.histogram(mt2.call_rate, range=(0, 1.0), bins=100, title='Variant Call Rate Histogram', legend='Call Rate')\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise: GQ vs DP\n",
    "\n",
    "In this exercise, you'll use Hail to investigate a strange property of sequencing datasets.\n",
    "\n",
    "The `DP` field is the sequencing depth (the number of reads).\n",
    "\n",
    "Let's first plot a histogram of `DP`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "p = hl.plot.histogram(mt.DP, range=(0, 40), bins=40, title='DP Histogram', legend='DP')\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now, let's do the same thing for GQ.\n",
    "\n",
    "The `GQ` field is the phred-scaled \"genotype quality\". The formula to convert to a linear-scale confidence (0 to 1) is  `10 ** -(mt.GQ / 10)`. GQ is truncated to lie between 0 and 99.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "p = hl.plot.histogram(mt.GQ, range=(0, 100), bins=100, title='GQ Histogram', legend='GQ')\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whoa!  That's a strange distribution!  There's a big spike at 100. The rest of the values have roughly the same shape as the DP distribution, but form a [Dimetrodon](https://en.wikipedia.org/wiki/Dimetrodon).  Use Hail to figure out what's going on!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
